########### MANUSCRIPT: The Role of Education and Attitudes in Cooking Fuel Choice: Evidence from two states in India
########### JOURNAL: ENERGY FOR SUSTAINABLE DEVELOPMENT
########### AUTHORS: CARLOS F. GOULD/1 AND JOHANNES URPELAINEN/2
########### AFFILIATIONS: 1/COLUMBIA UNIVERSITY MAILMAN SCHOOL OF PUBLIC HEALTH AND 2/JOHNS HOPKINS SCHOOL OF ADVANCED INTERNATIONAL STUDIES
########### PURPOSE: THIS IS CODE #5 THAT RUNS REGRESSIONS FOR THE MAIN TEXT


# for ordering covariates in the figures

covs_order <- c("Primary School", "Middle or High School", "Greater than High School",
                "Total monthly household expenditure\n (Logarithmized)",
                "Number of adults (>=18 years)", "Number of children (<18 years)",
                "Age of respondent (years)", 
                "Other backward caste", "Scheduled caste", "Scheduled tribe", 
                "Muslim", "Christian", "Religion other",
                "Rural (=1)")


# Figure 1 ----------------------------

# Coeffcient plot for logistic regressions showing the association between CWE education
# and LPG ownership in Kerala and Rajasthan. Regressions control for covariates and respective
# district fixed effects. Coeffcients are not exponentiated. Whiskers show 95% confidence intervals.
# Stars show statistical significance after Benjamini-Hochberg method for controlling for false
# discovery rate: * P<0.10, ** P<0.05, *** P<0.01.

# KERALA
reg2_K <- glm(Owns_LPG ~ 
                Education_CWE_primary + Education_CWE_middle_high + Education_CWE_greaterhigh +
                Exp_log + NumberAdults + NumberChildrenUnder18 +
                AgeRespondent + 
                Caste_OBC + Caste_ScheduledCaste + Caste_ScheduledTribe + Caste_DKCS +
                Religion_Muslim + Religion_Christian + Religion_Other + 
                Rural + 
                factor(District), data=kerala, family="binomial")

# RAJASTHAN
reg2_R <- glm(Owns_LPG ~ 
                Education_CWE_primary + Education_CWE_middle_high + Education_CWE_greaterhigh +
                Exp_log + NumberAdults + NumberChildrenUnder18 +
                AgeRespondent + 
                Caste_OBC + Caste_ScheduledCaste + Caste_ScheduledTribe + Caste_DKCS +
                Religion_Muslim + Religion_Christian + Religion_Other + 
                Rural + 
                factor(District), data=rajasthan, family="binomial")


DistrictFE <- c("factor(District)2", "factor(District)3", "factor(District)4", "factor(District)5", 
                "factor(District)6", "factor(District)7", "factor(District)8", 
                "factor(District)10", "factor(District)11", "factor(District)12", "factor(District)13", 
                "factor(District)14", "factor(District)15", "factor(District)16", "factor(District)17",
                "factor(District)18", "factor(District)19", "factor(District)20", "factor(District)21",
                "factor(District)22", "factor(District)23", "factor(District)24", "factor(District)25",
                "factor(District)26", "factor(District)27", "factor(District)28")

reg12_df <- rbind(tidy(reg2_K) %>% mutate(model="Kerala"),
                  tidy(reg2_R) %>% mutate(model="Rajasthan")) %>%
  filter(!(term %in% DistrictFE)) %>%
  filter(!(term %in% "Caste_DKCS")) %>%
  relabel_predictors(c(Education_CWE_primary="Primary School",
                       Education_CWE_middle_high="Middle or High School",
                       Education_CWE_greaterhigh="Greater than High School",
                       Exp_log="Total monthly household expenditure\n (Logarithmized)",
                       NumberAdults="Number of adults (>=18 years)",
                       NumberChildrenUnder18="Number of children (<18 years)",
                       AgeRespondent="Age of respondent (years)",
                       Caste_OBC="Other backward caste",
                       Caste_ScheduledCaste="Scheduled caste",
                       Caste_ScheduledTribe="Scheduled tribe",
                       Caste_DKCS="Caste DK/CS",
                       Religion_Muslim="Muslim",
                       Religion_Christian="Christian",
                       Religion_Other="Religion other",
                       Rural="Rural (=1)"))

reg12_df$p.value.bh <- p.adjust(reg12_df$p.value, method = "BH", n = length(reg12_df$p.value))
reg12_df$p.value.bh.star <- ifelse(reg12_df$p.value.bh<0.01, "***",
                                   ifelse(reg12_df$p.value.bh<0.05, "**",
                                          ifelse(reg12_df$p.value.bh<0.10, "*","")))

reg12_df[reg12_df=="Kerala"] <- "Kerala (With District FEs)"
reg12_df[reg12_df=="Rajasthan"] <- "Rajasthan (With District FEs)"


ggplot(reg12_df, aes(x=term, y=estimate,
                     ymin=estimate - 1.96*std.error, ymax = estimate + 1.96*std.error, label=p.value.bh.star,
                     color=model, shape=model)) + 
  geom_point(position = position_dodge(width = .5), size=3) +
  geom_errorbar(position = position_dodge(width = .5), size=1.5, width = 0) + 
  geom_text(aes(label=p.value.bh.star, x=term, y=((estimate+1.96*std.error)+0.2)), 
            position = position_dodge(width = .5), vjust=.70, size = 6, 
            show.legend = F) +
  geom_hline(yintercept=0, colour="grey60", linetype=2) + 
  scale_color_manual(values=c("black", "grey70")) +
  scale_shape_manual(values=c(16,15)) +
  scale_x_discrete(limits=rev(covs_order)) +
  theme_bw() +
  xlab("") + ylab("") +
  ggtitle("Association Between CWE Educational Achievement and LPG Ownership",
          subtitle = "Kerala (N=3929) and Rajasthan (N=6077)") +
  coord_flip(ylim = c(-2.5,5)) + 
  theme(legend.position = c(0.85,0.95),
        legend.title = element_blank(),
        legend.text = element_text(size=12),
        plot.title = element_text(size=18),
        plot.subtitle = element_text(size=16),
        axis.text = element_text(size=15),
        strip.text.x = element_text(size=16),
        strip.text.y = element_text(size=14)) 



# Figure 2 ----------------------------------

# Coeffcient plot for logistic regressions showing the association between the highest education
# level achieved by a female in the household and LPG ownership, controlling for covariates.
# Models 1 and 2 also include the caste category "Don't know / Can't say" (not shown on
# plot) and model 2 accounts for district fixed effects (not shown on plot). Coeffcients are not
# exponentiated. Whiskers show 95% confidence intervals. Stars show statistical significance after
# Benjamini-Hochberg method for controlling for false discovery rate: * P<0.10, ** P<0.05, ***
#   P<0.01.



# KERALA
reg1_f_k <- glm(Owns_LPG ~ 
                  Education_female_primary + Education_female_middle_high + Education_female_greaterhigh +
                  Exp_log + NumberAdults + NumberChildrenUnder18 +
                  AgeRespondent + 
                  Caste_OBC + Caste_ScheduledCaste + Caste_ScheduledTribe + Caste_DKCS +
                  Religion_Muslim + Religion_Christian + Religion_Other + 
                  Rural, data=kerala, family="binomial")

reg2_f_k <- glm(Owns_LPG ~ 
                  Education_female_primary + Education_female_middle_high + Education_female_greaterhigh +
                  Exp_log + NumberAdults + NumberChildrenUnder18 +
                  AgeRespondent + 
                  Caste_OBC + Caste_ScheduledCaste + Caste_ScheduledTribe + Caste_DKCS +
                  Religion_Muslim + Religion_Christian + Religion_Other + 
                  Rural + 
                  factor(District), data=kerala, family="binomial")

# RAJASTHAN
reg1_f_r <- glm(Owns_LPG ~ 
                  Education_4_FemaleSpouse_primary + Education_4_FemaleSpouse_middle_high + Education_4_FemaleSpouse_greaterhigh +
                  Exp_log + NumberAdults + NumberChildrenUnder18 +
                  AgeRespondent + 
                  Caste_OBC + Caste_ScheduledCaste + Caste_ScheduledTribe + Caste_DKCS +
                  Religion_Muslim + Religion_Christian + Religion_Other + 
                  Rural, data=rajasthan, family="binomial")


reg2_f_r <- glm(Owns_LPG ~ 
                  Education_4_FemaleSpouse_primary + Education_4_FemaleSpouse_middle_high + Education_4_FemaleSpouse_greaterhigh +
                  Exp_log + NumberAdults + NumberChildrenUnder18 +
                  AgeRespondent + 
                  Caste_OBC + Caste_ScheduledCaste + Caste_ScheduledTribe + Caste_DKCS +
                  Religion_Muslim + Religion_Christian + Religion_Other + 
                  Rural + 
                  factor(District), data=rajasthan, family="binomial")


DistrictFE <- c("factor(District)2", "factor(District)3", "factor(District)4", "factor(District)5", 
                "factor(District)6", "factor(District)7", "factor(District)8", 
                "factor(District)10", "factor(District)11", "factor(District)12", "factor(District)13", 
                "factor(District)14", "factor(District)15", "factor(District)16", "factor(District)17",
                "factor(District)18", "factor(District)19", "factor(District)20", "factor(District)21",
                "factor(District)22", "factor(District)23", "factor(District)24", "factor(District)25",
                "factor(District)26", "factor(District)27", "factor(District)28")

reg12_df <- rbind(tidy(reg1_f_k) %>% mutate(model="Model 1, No District FE",
                                            state = "Kerala"),
                  tidy(reg2_f_k) %>% mutate(model="Model 2, District FE",
                                            state = "Kerala"),
                  tidy(reg1_f_k) %>% mutate(model="Model 1, No District FE",
                                            state = "Rajasthan"),
                  tidy(reg2_f_r) %>% mutate(model="Model 2, District FE",
                                            state = "Rajasthan")) %>%
  filter(!(term %in% DistrictFE)) %>%
  filter(!(term %in% "Caste_DKCS")) %>%
  relabel_predictors(c(Education_female_primary="Primary School",
                       Education_female_middle_high="Middle or High School",
                       Education_female_greaterhigh="Greater than High School",
                       Education_4_FemaleSpouse_primary="Primary School",
                       Education_4_FemaleSpouse_middle_high="Middle or High School",
                       Education_4_FemaleSpouse_greaterhigh="Greater than High School",
                       Exp_log="Total monthly household expenditure\n (Logarithmized)",
                       NumberAdults="Number of adults (>=18 years)",
                       NumberChildrenUnder18="Number of children (<18 years)",
                       AgeRespondent="Age of respondent (years)",
                       Caste_OBC="Other backward caste",
                       Caste_ScheduledCaste="Scheduled caste",
                       Caste_ScheduledTribe="Scheduled tribe",
                       Caste_DKCS="Caste DK/CS",
                       Religion_Muslim="Muslim",
                       Religion_Christian="Christian",
                       Religion_Other="Religion other",
                       Rural="Rural (=1)"))

reg12_df$p.value.bh <- p.adjust(reg12_df$p.value, method = "BH", n = length(reg12_df$p.value))
reg12_df$p.value.bh.star <- ifelse(reg12_df$p.value.bh<0.01, "***",
                                   ifelse(reg12_df$p.value.bh<0.05, "**",
                                          ifelse(reg12_df$p.value.bh<0.10, "*","")))



ggplot(reg12_df, aes(x=term, y=estimate,
                     ymin=estimate - 1.96*std.error, ymax = estimate + 1.96*std.error, label=p.value.bh.star,
                     color=model, shape=model)) + 
  geom_point(position = position_dodge(width = .5), size=3) +
  geom_errorbar(position = position_dodge(width = .5), size=1.5, width = 0) + 
  geom_text(aes(label=p.value.bh.star, x=term, y=((estimate+1.96*std.error)+0.2)), 
            position = position_dodge(width = .5), vjust=.70, size = 6, 
            show.legend = F) +
  geom_hline(yintercept=0, colour="grey60", linetype=2) + 
  scale_color_manual(values=c("black", "grey70")) +
  scale_shape_manual(values=c(16,15)) +
  scale_x_discrete(limits=rev(covs_order)) +
  theme_bw() +
  xlab("") + ylab("") +
  ggtitle("Female Educational Achievement and LPG Ownership") +
  coord_flip(ylim = c(-2.5,5)) + 
  theme(legend.position = c(-0.22,1.04),
        legend.title = element_blank(),
        legend.text = element_text(size=16),
        plot.title = element_text(size=18),
        plot.subtitle = element_text(size=16),
        axis.text = element_text(size=15),
        strip.text.x = element_text(size=16),
        strip.text.y = element_text(size=14)) +
  facet_grid(.~state)




# Figure 4 ----------------------------------

# Coefficient plot for logistic regressions showing the association between perception indices
# and LPG ownership, controlling for covariates. Model 2 accounts for district fixed effects (not shown on plot). 
# Coefficients are not exponentiated. Whiskers show 95% confidence intervals. Stars show
# statistical significance after Benjamini-Hochberg method for controlling for false discovery rate: *
#   P<0.10, ** P<0.05, *** P<0.01.


k_preg1a <- glm(Owns_LPG ~ Perception_Index_FirewoodValuable_Norm +
                Exp_log + NumberAdults + NumberChildrenUnder18 +
                AgeRespondent + 
                Caste_OBC + Caste_ScheduledCaste + Caste_ScheduledTribe + Caste_DKCS +
                Religion_Muslim + Religion_Christian + Religion_Other + 
                Rural, data=kerala, family = "binomial")

k_preg1 <- glm(Owns_LPG ~ Perception_Index_FirewoodValuable_Norm +
               Exp_log + NumberAdults + NumberChildrenUnder18 +
               AgeRespondent + 
               Caste_OBC + Caste_ScheduledCaste + Caste_ScheduledTribe + Caste_DKCS +
               Religion_Muslim + Religion_Christian + Religion_Other + 
               Rural + 
               factor(District), data=kerala, family = "binomial")


k_preg2a <- glm(Owns_LPG ~ Perception_Index_ChulhaDifficult_Norm +
                Exp_log + NumberAdults + NumberChildrenUnder18 +
                AgeRespondent + 
                Caste_OBC + Caste_ScheduledCaste + Caste_ScheduledTribe + Caste_DKCS +
                Religion_Muslim + Religion_Christian + Religion_Other + 
                Rural, data=kerala, family = "binomial")

k_preg2 <- glm(Owns_LPG ~ Perception_Index_ChulhaDifficult_Norm +
               Exp_log + NumberAdults + NumberChildrenUnder18 +
               AgeRespondent + 
               Caste_OBC + Caste_ScheduledCaste + Caste_ScheduledTribe + Caste_DKCS +
               Religion_Muslim + Religion_Christian + Religion_Other + 
               Rural + 
               factor(District), data=kerala, family = "binomial")


k_preg3a <- glm(Owns_LPG ~ Perception_Index_FirewoodUnavailableDifficult_Norm +
                Exp_log + NumberAdults + NumberChildrenUnder18 +
                AgeRespondent + 
                Caste_OBC + Caste_ScheduledCaste + Caste_ScheduledTribe + Caste_DKCS +
                Religion_Muslim + Religion_Christian + Religion_Other + 
                Rural, data=kerala, family = "binomial")

k_preg3 <- glm(Owns_LPG ~ Perception_Index_FirewoodUnavailableDifficult_Norm +
               Exp_log + NumberAdults + NumberChildrenUnder18 +
               AgeRespondent + 
               Caste_OBC + Caste_ScheduledCaste + Caste_ScheduledTribe + Caste_DKCS +
               Religion_Muslim + Religion_Christian + Religion_Other + 
               Rural + 
               factor(District), data=kerala, family = "binomial")

k_preg4a <- glm(Owns_LPG ~ Perception_Index_CookingHealth_Norm +
                Exp_log + NumberAdults + NumberChildrenUnder18 +
                AgeRespondent + 
                Caste_OBC + Caste_ScheduledCaste + Caste_ScheduledTribe + Caste_DKCS +
                Religion_Muslim + Religion_Christian + Religion_Other + 
                Rural, data=kerala, family = "binomial")

k_preg4 <- glm(Owns_LPG ~ Perception_Index_CookingHealth_Norm +
               Exp_log + NumberAdults + NumberChildrenUnder18 +
               AgeRespondent + 
               Caste_OBC + Caste_ScheduledCaste + Caste_ScheduledTribe + Caste_DKCS +
               Religion_Muslim + Religion_Christian + Religion_Other + 
               Rural + 
               factor(District), data=kerala, family = "binomial")


k_preg5a <- glm(Owns_LPG ~ Perception_Index_GeneralHealth_Norm +
                Exp_log + NumberAdults + NumberChildrenUnder18 +
                AgeRespondent + 
                Caste_OBC + Caste_ScheduledCaste + Caste_ScheduledTribe + Caste_DKCS +
                Religion_Muslim + Religion_Christian + Religion_Other + 
                Rural + 
                factor(District), data=kerala, family = "binomial")

k_preg5 <- glm(Owns_LPG ~ Perception_Index_GeneralHealth_Norm +
               Exp_log + NumberAdults + NumberChildrenUnder18 +
               AgeRespondent + 
               Caste_OBC + Caste_ScheduledCaste + Caste_ScheduledTribe + Caste_DKCS +
               Religion_Muslim + Religion_Christian + Religion_Other + 
               Rural, data=kerala, family = "binomial")


k_preg6a <- glm(Owns_LPG ~ Perception_Index_LPGGoodAffordable_Norm +
                Exp_log + NumberAdults + NumberChildrenUnder18 +
                AgeRespondent + 
                Caste_OBC + Caste_ScheduledCaste + Caste_ScheduledTribe + Caste_DKCS +
                Religion_Muslim + Religion_Christian + Religion_Other + 
                Rural, data=kerala, family = "binomial")

k_preg6 <- glm(Owns_LPG ~ Perception_Index_LPGGoodAffordable_Norm +
               Exp_log + NumberAdults + NumberChildrenUnder18 +
               AgeRespondent + 
               Caste_OBC + Caste_ScheduledCaste + Caste_ScheduledTribe + Caste_DKCS +
               Religion_Muslim + Religion_Christian + Religion_Other + 
               Rural + 
               factor(District), data=kerala, family = "binomial")


k_preg_df1 <- data.frame(matrix(NA,1,5))
colnames(k_preg_df1) <- c("term", "estimate", "std.error", "statistic", "p.value")
k_preg_df <- rbind(k_preg_df1,
                 tidy(k_preg1a)[2,],
                 tidy(k_preg1)[2,],
                 tidy(k_preg2a)[2,],
                 tidy(k_preg2)[2,],
                 tidy(k_preg3a)[2,],
                 tidy(k_preg3)[2,],
                 tidy(k_preg4a)[2,],
                 tidy(k_preg4)[2,],
                 tidy(k_preg5a)[2,],
                 tidy(k_preg5)[2,],
                 tidy(k_preg6a)[2,],
                 tidy(k_preg6)[2,])
k_preg_df <- k_preg_df[-1,]
k_preg_df$model <- c("Model 1, No District FE", "Model 2, With FE")

k_preg_df[k_preg_df=="Perception_Index_LPGGoodAffordable_Norm"] <- "LPG is good and affordable"
k_preg_df[k_preg_df=="Perception_Index_FirewoodValuable_Norm"] <- "Chulha and firewood are valuable"
k_preg_df[k_preg_df=="Perception_Index_FirewoodUnavailableDifficult_Norm"] <- "Woodfuels are not available and inconvenient to use"
k_preg_df[k_preg_df=="Perception_Index_CookingHealth_Norm"] <- "Health attitudes in cooking"
k_preg_df[k_preg_df=="Perception_Index_GeneralHealth_Norm"] <- "General health attitudes are progressive"
k_preg_df[k_preg_df=="Perception_Index_ChulhaDifficult_Norm"] <- "Using the chulha is difficult"

k_preg_df$index <- k_preg_df$term
k_preg_df$p.value.bh <- p.adjust(k_preg_df$p.value, method = "BH", n = length(k_preg_df$p.value))
k_preg_df$p.value.bh.star <- ifelse(k_preg_df$p.value.bh<0.01, "***",
                                  ifelse(k_preg_df$p.value.bh<0.05, "**",
                                         ifelse(k_preg_df$p.value.bh<0.10, "*","")))
k_preg_df["State"] <- "Kerala"


# RAJASTHAN

r_preg1a <- glm(Owns_LPG ~ Perception_Index_LPGGood +
                Exp_log + NumberAdults + NumberChildrenUnder18 +
                AgeRespondent + 
                Caste_OBC + Caste_ScheduledCaste + Caste_ScheduledTribe + Caste_DKCS +
                Religion_Muslim + Religion_Christian + Religion_Other + 
                Rural, data=rajasthan, family = "binomial")

r_preg1 <- glm(Owns_LPG ~ Perception_Index_LPGGood +
               Exp_log + NumberAdults + NumberChildrenUnder18 +
               AgeRespondent + 
               Caste_OBC + Caste_ScheduledCaste + Caste_ScheduledTribe + Caste_DKCS +
               Religion_Muslim + Religion_Christian + Religion_Other + 
               Rural + 
               factor(District), data=rajasthan, family = "binomial")


r_preg2a <- glm(Owns_LPG ~ Perception_Index_LPGaffordable +
                Exp_log + NumberAdults + NumberChildrenUnder18 +
                AgeRespondent + 
                Caste_OBC + Caste_ScheduledCaste + Caste_ScheduledTribe + Caste_DKCS +
                Religion_Muslim + Religion_Christian + Religion_Other + 
                Rural, data=rajasthan, family = "binomial")

r_preg2 <- glm(Owns_LPG ~ Perception_Index_LPGaffordable +
               Exp_log + NumberAdults + NumberChildrenUnder18 +
               AgeRespondent + 
               Caste_OBC + Caste_ScheduledCaste + Caste_ScheduledTribe + Caste_DKCS +
               Religion_Muslim + Religion_Christian + Religion_Other + 
               Rural + 
               factor(District), data=rajasthan, family = "binomial")


r_preg3a <- glm(Owns_LPG ~ Perception_Index_CookingHealth +
                Exp_log + NumberAdults + NumberChildrenUnder18 +
                AgeRespondent + 
                Caste_OBC + Caste_ScheduledCaste + Caste_ScheduledTribe + Caste_DKCS +
                Religion_Muslim + Religion_Christian + Religion_Other + 
                Rural, data=rajasthan, family = "binomial")

r_preg3 <- glm(Owns_LPG ~ Perception_Index_CookingHealth +
               Exp_log + NumberAdults + NumberChildrenUnder18 +
               AgeRespondent + 
               Caste_OBC + Caste_ScheduledCaste + Caste_ScheduledTribe + Caste_DKCS +
               Religion_Muslim + Religion_Christian + Religion_Other + 
               Rural + 
               factor(District), data=rajasthan, family = "binomial")

r_preg4a <- glm(Owns_LPG ~ Perception_Index_FirewoodValuable +
                Exp_log + NumberAdults + NumberChildrenUnder18 +
                AgeRespondent + 
                Caste_OBC + Caste_ScheduledCaste + Caste_ScheduledTribe + Caste_DKCS +
                Religion_Muslim + Religion_Christian + Religion_Other + 
                Rural, data=rajasthan, family = "binomial")

r_preg4 <- glm(Owns_LPG ~ Perception_Index_FirewoodValuable +
               Exp_log + NumberAdults + NumberChildrenUnder18 +
               AgeRespondent + 
               Caste_OBC + Caste_ScheduledCaste + Caste_ScheduledTribe + Caste_DKCS +
               Religion_Muslim + Religion_Christian + Religion_Other + 
               Rural + 
               factor(District), data=rajasthan, family = "binomial")


r_preg5a <- glm(Owns_LPG ~ Perception_Index_FirewoodUnavailableDifficult +
                Exp_log + NumberAdults + NumberChildrenUnder18 +
                AgeRespondent + 
                Caste_OBC + Caste_ScheduledCaste + Caste_ScheduledTribe + Caste_DKCS +
                Religion_Muslim + Religion_Christian + Religion_Other + 
                Rural + 
                factor(District), data=rajasthan, family = "binomial")

r_preg5 <- glm(Owns_LPG ~ Perception_Index_FirewoodUnavailableDifficult +
               Exp_log + NumberAdults + NumberChildrenUnder18 +
               AgeRespondent + 
               Caste_OBC + Caste_ScheduledCaste + Caste_ScheduledTribe + Caste_DKCS +
               Religion_Muslim + Religion_Christian + Religion_Other + 
               Rural, data=rajasthan, family = "binomial")



r_preg_df1 <- data.frame(matrix(NA,1,5))
colnames(r_preg_df1) <- c("term", "estimate", "std.error", "statistic", "p.value")
r_preg_df <- rbind(r_preg_df1,
                 tidy(r_preg1a)[2,],
                 tidy(r_preg1)[2,],
                 tidy(r_preg2a)[2,],
                 tidy(r_preg2)[2,],
                 tidy(r_preg3a)[2,],
                 tidy(r_preg3)[2,],
                 tidy(r_preg4a)[2,],
                 tidy(r_preg4)[2,],
                 tidy(r_preg5a)[2,],
                 tidy(r_preg5)[2,])
r_preg_df <- r_preg_df[-1,]
r_preg_df$model <- c("Model 1, No District FE", "Model 2, With FE")

r_preg_df$index <- r_preg_df$term
r_preg_df$p.value.bh <- p.adjust(r_preg_df$p.value, method = "BH", n = length(r_preg_df$p.value))
r_preg_df$p.value.bh.star <- ifelse(r_preg_df$p.value.bh<0.01, "***",
                                  ifelse(r_preg_df$p.value.bh<0.05, "**",
                                         ifelse(r_preg_df$p.value.bh<0.10, "*","")))
r_preg_df["State"] <- "Rajasthan"


r_preg_df[r_preg_df=="Perception_Index_FirewoodUnavailableDifficult"] <- "Woodfuels are not available, inconvenient"
r_preg_df[r_preg_df=="Perception_Index_LPGGood"] <- "LPG is good quality cooking and beneficial"
r_preg_df[r_preg_df=="Perception_Index_LPGaffordable"] <- "LPG is desirable, but expensive"
r_preg_df[r_preg_df=="Perception_Index_CookingHealth"] <- "Cooking fuel choices affect health"
r_preg_df[r_preg_df=="Perception_Index_FirewoodValuable"] <- "Chulha and woodfuel are valuable"

combine_indices_lpg <- k_preg_df %>%
  bind_rows(r_preg_df)

combine_indices_lpg <- combine_indices_lpg[-1,]

combine_indices_lpg[combine_indices_lpg=="Chulha and firewood are valuable"] <- "Chulha and woodfuel are valuable"
combine_indices_lpg[combine_indices_lpg=="Health attitudes in cooking"] <- "Cooking fuel choices affect health"


ggplot(combine_indices_lpg, aes(x=index, y=estimate,
                                ymin=estimate - 1.96*std.error, ymax = estimate + 1.96*std.error, label=p.value.bh.star,
                                shape=model, color=model)) + 
  geom_point(position = position_dodge(width = .5), size=3) +
  geom_errorbar(position = position_dodge(width = .5), size=1.5, width = 0) + 
  geom_text(aes(label=p.value.bh.star, x=index, y=((estimate - 1.96*std.error)-0.2)), position = position_dodge(width = .5), vjust=0.68, size = 6) +
  geom_hline(yintercept=0, colour="grey60", linetype=2) + 
  theme_bw() +
  xlab("") + ylab("") +
  scale_shape_manual(values = c(16, 15), guide = FALSE) +
  scale_linetype_manual(values = c(1,5), guide = FALSE) +
  scale_color_manual(values = c("black", "grey60"), guide = guide_legend(reverse = TRUE)) +
  ggtitle("Perception Indices and LPG Ownership") +
  coord_flip() + 
  guides(colour = guide_legend(reverse=T)) +
  theme(legend.position = c(-0.35,1.04),
        legend.title = element_blank(),
        legend.text = element_text(size=12),
        legend.background = element_blank(),
        legend.box.background = element_rect(colour = "black"),
        plot.title = element_text(size=16),
        axis.text = element_text(size=15),
        strip.text.x = element_text(size=14),
        strip.text.y = element_text(size=14)) +
  facet_grid(.~State)





# Figure 5 -------------------------------

### Figure 5 Kerala -----

# Education predicting perception "firewood is valuable"
k_ed_Perception1 <- lm(Perception_Index_FirewoodValuable_Norm ~ 
                    Education_CWE_primary + Education_CWE_middle_high + Education_CWE_greaterhigh +
                    Exp_log + NumberAdults + NumberChildrenUnder18 +
                    AgeRespondent + 
                    Caste_OBC + Caste_ScheduledCaste + Caste_ScheduledTribe + Caste_DKCS +
                    Religion_Muslim + Religion_Christian + Religion_Other + 
                    Rural + factor(District), data=kerala)

# Education predicting perception "Chulha is difficult to use"
k_ed_Perception2 <- lm(Perception_Index_ChulhaDifficult_Norm ~ 
                    Education_CWE_primary + Education_CWE_middle_high + Education_CWE_greaterhigh +
                    Exp_log + NumberAdults + NumberChildrenUnder18 +
                    AgeRespondent + 
                    Caste_OBC + Caste_ScheduledCaste + Caste_ScheduledTribe + Caste_DKCS +
                    Religion_Muslim + Religion_Christian + Religion_Other + 
                    Rural + factor(District), data=kerala)

# Education predicting perception "Firewood is unavailable and/or difficult to acquire"

k_ed_Perception3 <- lm(Perception_Index_FirewoodUnavailableDifficult_Norm ~ 
                    Education_CWE_primary + Education_CWE_middle_high + Education_CWE_greaterhigh +
                    Exp_log + NumberAdults + NumberChildrenUnder18 +
                    AgeRespondent + 
                    Caste_OBC + Caste_ScheduledCaste + Caste_ScheduledTribe + Caste_DKCS +
                    Religion_Muslim + Religion_Christian + Religion_Other + 
                    Rural + factor(District), data=kerala)

# Education predicting perception "Cooking fuel has a relationship with health"

k_ed_Perception4 <- lm(Perception_Index_CookingHealth_Norm ~ 
                    Education_CWE_primary + Education_CWE_middle_high + Education_CWE_greaterhigh +
                    Exp_log + NumberAdults + NumberChildrenUnder18 +
                    AgeRespondent + 
                    Caste_OBC + Caste_ScheduledCaste + Caste_ScheduledTribe + Caste_DKCS +
                    Religion_Muslim + Religion_Christian + Religion_Other + 
                    Rural + factor(District), data=kerala)

# Education predicting perception "General progressive health norms"

k_ed_Perception5 <- lm(Perception_Index_GeneralHealth_Norm ~ 
                    Education_CWE_primary + Education_CWE_middle_high + Education_CWE_greaterhigh +
                    Exp_log + NumberAdults + NumberChildrenUnder18 +
                    AgeRespondent + 
                    Caste_OBC + Caste_ScheduledCaste + Caste_ScheduledTribe + Caste_DKCS +
                    Religion_Muslim + Religion_Christian + Religion_Other + 
                    Rural + factor(District), data=kerala)

# Education predicting perception "LPG is a good fuel and is affordable"

k_ed_Perception6 <- lm(Perception_Index_LPGGoodAffordable_Norm ~ 
                    Education_CWE_primary + Education_CWE_middle_high + Education_CWE_greaterhigh +
                    Exp_log + NumberAdults + NumberChildrenUnder18 +
                    AgeRespondent + 
                    Caste_OBC + Caste_ScheduledCaste + Caste_ScheduledTribe + Caste_DKCS +
                    Religion_Muslim + Religion_Christian + Religion_Other + 
                    Rural + factor(District), data=kerala)

# Female education

# Education predicting Perception "firewood is valuable"
k_f_ed_Perception1 <- lm(Perception_Index_FirewoodValuable_Norm ~ 
                    Education_female_primary + Education_female_middle_high + Education_female_greaterhigh +
                    Exp_log + NumberAdults + NumberChildrenUnder18 +
                    AgeRespondent + 
                    Caste_OBC + Caste_ScheduledCaste + Caste_ScheduledTribe + Caste_DKCS +
                    Religion_Muslim + Religion_Christian + Religion_Other + 
                    Rural + factor(District), data=kerala)

# Education predicting Perception "Chulha is difficult to use"
k_f_ed_Perception2 <- lm(Perception_Index_ChulhaDifficult_Norm ~ 
                    Education_female_primary + Education_female_middle_high + Education_female_greaterhigh +
                    Exp_log + NumberAdults + NumberChildrenUnder18 +
                    AgeRespondent + 
                    Caste_OBC + Caste_ScheduledCaste + Caste_ScheduledTribe + Caste_DKCS +
                    Religion_Muslim + Religion_Christian + Religion_Other + 
                    Rural + factor(District), data=kerala)

# Education predicting Perception "Firewood is unavailable and/or difficult to acquire"

k_f_ed_Perception3 <- lm(Perception_Index_FirewoodUnavailableDifficult_Norm ~ 
                    Education_female_primary + Education_female_middle_high + Education_female_greaterhigh +
                    Exp_log + NumberAdults + NumberChildrenUnder18 +
                    AgeRespondent + 
                    Caste_OBC + Caste_ScheduledCaste + Caste_ScheduledTribe + Caste_DKCS +
                    Religion_Muslim + Religion_Christian + Religion_Other + 
                    Rural + factor(District), data=kerala)

# Education predicting Perception "Cooking fuel has a relationship with health"

k_f_ed_Perception4 <- lm(Perception_Index_CookingHealth_Norm ~ 
                    Education_female_primary + Education_female_middle_high + Education_female_greaterhigh +
                    Exp_log + NumberAdults + NumberChildrenUnder18 +
                    AgeRespondent + 
                    Caste_OBC + Caste_ScheduledCaste + Caste_ScheduledTribe + Caste_DKCS +
                    Religion_Muslim + Religion_Christian + Religion_Other + 
                    Rural + factor(District), data=kerala)

# Education predicting Perception "General progressive health norms"

k_f_ed_Perception5 <- lm(Perception_Index_GeneralHealth_Norm ~ 
                    Education_female_primary + Education_female_middle_high + Education_female_greaterhigh +
                    Exp_log + NumberAdults + NumberChildrenUnder18 +
                    AgeRespondent + 
                    Caste_OBC + Caste_ScheduledCaste + Caste_ScheduledTribe + Caste_DKCS +
                    Religion_Muslim + Religion_Christian + Religion_Other + 
                    Rural + factor(District), data=kerala)

# Education predicting Perception "LPG is a good fuel and is affordable"

k_f_ed_Perception6 <- lm(Perception_Index_LPGGoodAffordable_Norm ~ 
                    Education_female_primary + Education_female_middle_high + Education_female_greaterhigh +
                    Exp_log + NumberAdults + NumberChildrenUnder18 +
                    AgeRespondent + 
                    Caste_OBC + Caste_ScheduledCaste + Caste_ScheduledTribe + Caste_DKCS +
                    Religion_Muslim + Religion_Christian + Religion_Other + 
                    Rural + factor(District), data=kerala)


k_ed_perception_df1 <- data.frame(matrix(NA,1,5))
colnames(k_ed_perception_df1) <- c("term", "estimate", "std.error", "statistic", "p.value")
k_ed_perception_df <- rbind(k_ed_perception_df1,
                            tidy(k_ed_Perception1)[2:4,],
                            tidy(k_f_ed_Perception1)[2:4,],
                            tidy(k_ed_Perception2)[2:4,],
                            tidy(k_f_ed_Perception2)[2:4,],
                            tidy(k_ed_Perception3)[2:4,],
                            tidy(k_f_ed_Perception3)[2:4,],
                            tidy(k_ed_Perception4)[2:4,],
                            tidy(k_f_ed_Perception4)[2:4,],
                            tidy(k_ed_Perception5)[2:4,],
                            tidy(k_f_ed_Perception5)[2:4,],
                            tidy(k_ed_Perception6)[2:4,],
                            tidy(k_f_ed_Perception6)[2:4,])
k_ed_perception_df <- k_ed_perception_df[-1,]
k_ed_perception_df$model <- c(rep("CWE",3), rep("Female",3))

k_ed_perception_df$index <- c(rep("Perception_Index_FirewoodValuable",6),
                              rep("Perception_Index_ChulhaDifficult",6),
                              rep("Perception_Index_FirewoodUnavailableDifficult",6),
                              rep("Perception_Index_CookingHealth",6),
                              rep("Perception_Index_GeneralHealth",6),
                              rep("Perception_Index_LPGGoodAffordable",6))

k_ed_perception_df$p.value.bh <- p.adjust(k_ed_perception_df$p.value, method = "BH", n = length(k_ed_perception_df$p.value))
k_ed_perception_df$p.value.bh.star <- ifelse(k_ed_perception_df$p.value.bh<0.01, "***",
                                             ifelse(k_ed_perception_df$p.value.bh<0.05, "**",
                                                    ifelse(k_ed_perception_df$p.value.bh<0.10, "*","")))

k_ed_perception_df[k_ed_perception_df=="Perception_Index_FirewoodValuable"] <- "Firewood is valuable"
k_ed_perception_df[k_ed_perception_df=="Perception_Index_ChulhaDifficult"] <- "Chulha is difficult to use"
k_ed_perception_df[k_ed_perception_df=="Perception_Index_CookingHealth"] <- "Cooking fuel type is related to health"
k_ed_perception_df[k_ed_perception_df=="Perception_Index_FirewoodUnavailableDifficult"] <- "Firewood is unavailable, inconvenient"
k_ed_perception_df[k_ed_perception_df=="Perception_Index_GeneralHealth"] <- "Generally progressive health norms"
k_ed_perception_df[k_ed_perception_df=="Perception_Index_LPGGoodAffordable"] <- "LPG is a good cooking fuel and affordable"

  
### Figure 5 Rajasthan -----


r_ed_Perception1 <- lm(Perception_Index_FirewoodUnavailableDifficult ~ 
                      Education_CWE_primary + Education_CWE_middle_high + Education_CWE_greaterhigh +
                      Exp_log + NumberAdults + NumberChildrenUnder18 +
                      AgeRespondent + 
                      Caste_OBC + Caste_ScheduledCaste + Caste_ScheduledTribe + Caste_DKCS +
                      Religion_Muslim + Religion_Christian + Religion_Other + 
                      Rural + factor(District), data=rajasthan)

r_ed_Perception2 <- lm(Perception_Index_LPGGood ~ 
                      Education_CWE_primary + Education_CWE_middle_high + Education_CWE_greaterhigh +
                      Exp_log + NumberAdults + NumberChildrenUnder18 +
                      AgeRespondent + 
                      Caste_OBC + Caste_ScheduledCaste + Caste_ScheduledTribe + Caste_DKCS +
                      Religion_Muslim + Religion_Christian + Religion_Other + 
                      Rural + factor(District), data=rajasthan)

r_ed_Perception3 <- lm(Perception_Index_LPGaffordable ~ 
                      Education_CWE_primary + Education_CWE_middle_high + Education_CWE_greaterhigh +
                      Exp_log + NumberAdults + NumberChildrenUnder18 +
                      AgeRespondent + 
                      Caste_OBC + Caste_ScheduledCaste + Caste_ScheduledTribe + Caste_DKCS +
                      Religion_Muslim + Religion_Christian + Religion_Other + 
                      Rural + factor(District), data=rajasthan)

r_ed_Perception4 <- lm(Perception_Index_CookingHealth ~ 
                      Education_CWE_primary + Education_CWE_middle_high + Education_CWE_greaterhigh +
                      Exp_log + NumberAdults + NumberChildrenUnder18 +
                      AgeRespondent + 
                      Caste_OBC + Caste_ScheduledCaste + Caste_ScheduledTribe + Caste_DKCS +
                      Religion_Muslim + Religion_Christian + Religion_Other + 
                      Rural + factor(District), data=rajasthan)

r_ed_Perception5 <- lm(Perception_Index_FirewoodValuable ~ 
                      Education_CWE_primary + Education_CWE_middle_high + Education_CWE_greaterhigh +
                      Exp_log + NumberAdults + NumberChildrenUnder18 +
                      AgeRespondent + 
                      Caste_OBC + Caste_ScheduledCaste + Caste_ScheduledTribe + Caste_DKCS +
                      Religion_Muslim + Religion_Christian + Religion_Other + 
                      Rural + factor(District), data=rajasthan)


# Rajasthan female spouse education

r_f_ed_Perception1 <- lm(Perception_Index_FirewoodUnavailableDifficult ~ 
                      Education_4_FemaleSpouse_primary + Education_4_FemaleSpouse_middle_high + Education_4_FemaleSpouse_greaterhigh +
                      Exp_log + NumberAdults + NumberChildrenUnder18 +
                      AgeRespondent + 
                      Caste_OBC + Caste_ScheduledCaste + Caste_ScheduledTribe + Caste_DKCS +
                      Religion_Muslim + Religion_Christian + Religion_Other + 
                      Rural + factor(District), data=rajasthan)

r_f_ed_Perception2 <- lm(Perception_Index_LPGGood ~ 
                      Education_4_FemaleSpouse_primary + Education_4_FemaleSpouse_middle_high + Education_4_FemaleSpouse_greaterhigh +
                      Exp_log + NumberAdults + NumberChildrenUnder18 +
                      AgeRespondent + 
                      Caste_OBC + Caste_ScheduledCaste + Caste_ScheduledTribe + Caste_DKCS +
                      Religion_Muslim + Religion_Christian + Religion_Other + 
                      Rural + factor(District), data=rajasthan)

r_f_ed_Perception3 <- lm(Perception_Index_LPGaffordable ~ 
                      Education_4_FemaleSpouse_primary + Education_4_FemaleSpouse_middle_high + Education_4_FemaleSpouse_greaterhigh +
                      Exp_log + NumberAdults + NumberChildrenUnder18 +
                      AgeRespondent + 
                      Caste_OBC + Caste_ScheduledCaste + Caste_ScheduledTribe + Caste_DKCS +
                      Religion_Muslim + Religion_Christian + Religion_Other + 
                      Rural + factor(District), data=rajasthan)

r_f_ed_Perception4 <- lm(Perception_Index_CookingHealth ~ 
                      Education_4_FemaleSpouse_primary + Education_4_FemaleSpouse_middle_high + Education_4_FemaleSpouse_greaterhigh +
                      Exp_log + NumberAdults + NumberChildrenUnder18 +
                      AgeRespondent + 
                      Caste_OBC + Caste_ScheduledCaste + Caste_ScheduledTribe + Caste_DKCS +
                      Religion_Muslim + Religion_Christian + Religion_Other + 
                      Rural + factor(District), data=rajasthan)

r_f_ed_Perception5 <- lm(Perception_Index_FirewoodValuable ~ 
                      Education_4_FemaleSpouse_primary + Education_4_FemaleSpouse_middle_high + Education_4_FemaleSpouse_greaterhigh +
                      Exp_log + NumberAdults + NumberChildrenUnder18 +
                      AgeRespondent + 
                      Caste_OBC + Caste_ScheduledCaste + Caste_ScheduledTribe + Caste_DKCS +
                      Religion_Muslim + Religion_Christian + Religion_Other + 
                      Rural + factor(District), data=rajasthan)


r_ed_perception_df1 <- data.frame(matrix(NA,1,5))
colnames(r_ed_perception_df1) <- c("term", "estimate", "std.error", "statistic", "p.value")
r_ed_perception_df <- rbind(r_ed_perception_df1,
                            tidy(r_ed_Perception1)[2:4,],
                            tidy(r_f_ed_Perception1)[2:4,],
                            tidy(r_ed_Perception2)[2:4,],
                            tidy(r_f_ed_Perception2)[2:4,],
                            tidy(r_ed_Perception3)[2:4,],
                            tidy(r_f_ed_Perception3)[2:4,],
                            tidy(r_ed_Perception4)[2:4,],
                            tidy(r_f_ed_Perception4)[2:4,],
                            tidy(r_ed_Perception5)[2:4,],
                            tidy(r_f_ed_Perception5)[2:4,])
r_ed_perception_df <- r_ed_perception_df[-1,]
r_ed_perception_df$model <- c(rep("CWE",3), rep("Female Spouse",3))

r_ed_perception_df$index <- c(rep("Perception_Index_LPGGood", 6),
                              rep("Perception_Index_LPGaffordable",6),
                              rep("Perception_Index_CookingHealth",6),
                              rep("Perception_Index_FirewoodValuable",6),
                              rep("Perception_Index_FirewoodUnavailableDifficult",6))

r_ed_perception_df$p.value.bh <- p.adjust(r_ed_perception_df$p.value, method = "BH", n = length(r_ed_perception_df$p.value))
r_ed_perception_df$p.value.bh.star <- ifelse(r_ed_perception_df$p.value.bh<0.01, "***",
                                             ifelse(r_ed_perception_df$p.value.bh<0.05, "**",
                                                    ifelse(r_ed_perception_df$p.value.bh<0.10, "*","")))

r_ed_perception_df[r_ed_perception_df=="Perception_Index_FirewoodUnavailableDifficult"] <- "Woodfuels are not available, inconvenient"
r_ed_perception_df[r_ed_perception_df=="Perception_Index_LPGGood"] <- "LPG is good quality cooking and beneficial"
r_ed_perception_df[r_ed_perception_df=="Perception_Index_LPGaffordable"] <- "LPG is desirable, but expensive"
r_ed_perception_df[r_ed_perception_df=="Perception_Index_CookingHealth"] <- "Cooking fuel choices affect health"
r_ed_perception_df[r_ed_perception_df=="Perception_Index_FirewoodValuable"] <- "Chulha and woodfuel are valuable"



# Figure 5: Making figure ------


k_ed_perception_df[k_ed_perception_df=="Education_CWE_primary"] <- "Primary School"
k_ed_perception_df[k_ed_perception_df=="Education_CWE_middle_high"] <- "Middle and High School"
k_ed_perception_df[k_ed_perception_df=="Education_CWE_greaterhigh"] <- "Greather than High School"

k_ed_perception_df[k_ed_perception_df=="Education_female_primary"] <- "Primary School"
k_ed_perception_df[k_ed_perception_df=="Education_female_middle_high"] <- "Middle and High School"
k_ed_perception_df[k_ed_perception_df=="Education_female_greaterhigh"] <- "Greather than High School"

k_index_ed_fig <-  ggplot(data=k_ed_perception_df, 
                          aes(x=factor(index), y=estimate,
                              ymin=estimate - 1.96*std.error, ymax = estimate + 1.96*std.error, 
                              label=p.value.bh.star,
                              shape=factor(term))) + 
  geom_point(position = position_dodge(width = .5), size=3) +
  geom_errorbar(position = position_dodge(width = .5), size=1.5, width = 0) + 
  geom_text(aes(label=p.value.bh.star, x = index, y = estimate - 2.25*std.error), 
            position = position_dodge(width = .5),
            vjust=0.65, size = 6) +
  geom_hline(yintercept=0, colour="grey60", linetype=2) + 
  theme_bw() +
  xlab("") + 
  ylab("") +
  ggtitle("Kerala: Perception Indices and Educational Attainment") +
  coord_flip() + 
  theme(legend.position = c(-0.35,1.04),
        legend.title = element_blank(),
        legend.text = element_text(size=12),
        legend.background = element_blank(),
        legend.box.background = element_rect(colour = "black"),
        plot.title = element_text(size=16),
        axis.text = element_text(size=15),
        strip.text.x = element_text(size=14),
        strip.text.y = element_text(size=14)) +
  facet_grid(.~model) 

k_index_ed_fig

r_ed_perception_df[r_ed_perception_df=="Education_CWE_primary"] <- "Primary School"
r_ed_perception_df[r_ed_perception_df=="Education_CWE_middle_high"] <- "Middle and High School"
r_ed_perception_df[r_ed_perception_df=="Education_CWE_greaterhigh"] <- "Greather than High School"

r_ed_perception_df[r_ed_perception_df=="Education_4_FemaleSpouse_primary"] <- "Primary School"
r_ed_perception_df[r_ed_perception_df=="Education_4_FemaleSpouse_middle_high"] <- "Middle and High School"
r_ed_perception_df[r_ed_perception_df=="Education_4_FemaleSpouse_greaterhigh"] <- "Greather than High School"



r_index_ed_fig <- ggplot(r_ed_perception_df, 
                         aes(x=factor(index), y=estimate,
                             ymin=estimate - 1.96*std.error, ymax = estimate + 1.96*std.error, label=p.value.bh.star,
                             shape=term)) + 
  geom_point(position = position_dodge(width = .5), size=3) +
  geom_errorbar(position = position_dodge(width = .5), size=1.5, width = 0) + 
  geom_text(aes(label=p.value.bh.star, x = index, y = estimate - 2.25*std.error), 
            position = position_dodge(width = .5),
            vjust=0.65, size = 6) +
  geom_hline(yintercept=0, colour="grey60", linetype=2) + 
  theme_bw() +
  xlab("") + ylab("") +
  ggtitle("Rajasthan: Perception Indices and Educational Attainment") +
  coord_flip() + 
  theme(legend.position = c(-0.35,1.04),
        legend.title = element_blank(),
        legend.text = element_text(size=12),
        legend.background = element_blank(),
        legend.box.background = element_rect(colour = "black"),
        plot.title = element_text(size=16),
        axis.text = element_text(size=15),
        strip.text.x = element_text(size=14),
        strip.text.y = element_text(size=14)) +
  facet_grid(.~model) 


plot_grid(
  
  k_index_ed_fig, 
  r_index_ed_fig,
  nrow=2
  
)

# Figure 6 --------------------------------------------

# Coefficient plot for linear regressions showing the associations between CWE and female
# education and LPG ownership, with and without including perception indices, in Kerala and Rajasthan
# households. Regressions also control for covariates and district fixed effects (not shown).
# Coefficients are not exponentiated. Whiskers show 95% confidence intervals. Stars show statistical
# significance after Benjamini-Hochberg method for controlling for false discovery rate: * P<0.10, **
#   P<0.05, *** P<0.01.

# Figure 6: Kerala CWE -----



mreg1a <- glm(Owns_LPG ~ 
                Education_CWE_primary + Education_CWE_middle_high + Education_CWE_greaterhigh +
                Exp_log + NumberAdults + NumberChildrenUnder18 +
                AgeRespondent + 
                Caste_OBC + Caste_ScheduledCaste + Caste_ScheduledTribe + Caste_DKCS +
                Religion_Muslim + Religion_Christian + Religion_Other + 
                Rural + factor(District), data=kerala, family = "binomial")

mreg1 <- glm(Owns_LPG ~ Perception_Index_FirewoodValuable_Norm +
               Education_CWE_primary + Education_CWE_middle_high + Education_CWE_greaterhigh +
               Exp_log + NumberAdults + NumberChildrenUnder18 +
               AgeRespondent + 
               Caste_OBC + Caste_ScheduledCaste + Caste_ScheduledTribe + Caste_DKCS +
               Religion_Muslim + Religion_Christian + Religion_Other + 
               Rural + 
               factor(District), data=kerala, family = "binomial")


mreg2a <- glm(Owns_LPG ~ 
                Education_CWE_primary + Education_CWE_middle_high + Education_CWE_greaterhigh +
                Exp_log + NumberAdults + NumberChildrenUnder18 +
                AgeRespondent + 
                Caste_OBC + Caste_ScheduledCaste + Caste_ScheduledTribe + Caste_DKCS +
                Religion_Muslim + Religion_Christian + Religion_Other + 
                Rural + factor(District), data=kerala, family = "binomial")

mreg2 <- glm(Owns_LPG ~ Perception_Index_ChulhaDifficult_Norm + 
               Education_CWE_primary + Education_CWE_middle_high + Education_CWE_greaterhigh +
               Exp_log + NumberAdults + NumberChildrenUnder18 +
               AgeRespondent + 
               Caste_OBC + Caste_ScheduledCaste + Caste_ScheduledTribe + Caste_DKCS +
               Religion_Muslim + Religion_Christian + Religion_Other + 
               Rural + 
               factor(District), data=kerala, family = "binomial")


mreg3a <- glm(Owns_LPG ~ 
                Education_CWE_primary + Education_CWE_middle_high + Education_CWE_greaterhigh +
                Exp_log + NumberAdults + NumberChildrenUnder18 +
                AgeRespondent + 
                Caste_OBC + Caste_ScheduledCaste + Caste_ScheduledTribe + Caste_DKCS +
                Religion_Muslim + Religion_Christian + Religion_Other + 
                Rural + factor(District), data=kerala, family = "binomial")

mreg3 <- glm(Owns_LPG ~ Perception_Index_FirewoodUnavailableDifficult_Norm +
               Education_CWE_primary + Education_CWE_middle_high + Education_CWE_greaterhigh +
               Exp_log + NumberAdults + NumberChildrenUnder18 +
               AgeRespondent + 
               Caste_OBC + Caste_ScheduledCaste + Caste_ScheduledTribe + Caste_DKCS +
               Religion_Muslim + Religion_Christian + Religion_Other + 
               Rural + 
               factor(District), data=kerala, family = "binomial")

mreg4a <- glm(Owns_LPG ~ 
                Education_CWE_primary + Education_CWE_middle_high + Education_CWE_greaterhigh +
                Exp_log + NumberAdults + NumberChildrenUnder18 +
                AgeRespondent + 
                Caste_OBC + Caste_ScheduledCaste + Caste_ScheduledTribe + Caste_DKCS +
                Religion_Muslim + Religion_Christian + Religion_Other + 
                Rural + factor(District), data=kerala, family = "binomial")

mreg4 <- glm(Owns_LPG ~ Perception_Index_CookingHealth_Norm +
               Education_CWE_primary + Education_CWE_middle_high + Education_CWE_greaterhigh +
               Exp_log + NumberAdults + NumberChildrenUnder18 +
               AgeRespondent + 
               Caste_OBC + Caste_ScheduledCaste + Caste_ScheduledTribe + Caste_DKCS +
               Religion_Muslim + Religion_Christian + Religion_Other + 
               Rural + 
               factor(District), data=kerala, family = "binomial")


mreg5a <- glm(Owns_LPG ~ 
                Education_CWE_primary + Education_CWE_middle_high + Education_CWE_greaterhigh +
                Exp_log + NumberAdults + NumberChildrenUnder18 +
                AgeRespondent + 
                Caste_OBC + Caste_ScheduledCaste + Caste_ScheduledTribe + Caste_DKCS +
                Religion_Muslim + Religion_Christian + Religion_Other + 
                Rural + 
                factor(District), data=kerala, family = "binomial")

mreg5 <- glm(Owns_LPG ~ Perception_Index_GeneralHealth_Norm +
               Education_CWE_primary + Education_CWE_middle_high + Education_CWE_greaterhigh +
               Exp_log + NumberAdults + NumberChildrenUnder18 +
               AgeRespondent + 
               Caste_OBC + Caste_ScheduledCaste + Caste_ScheduledTribe + Caste_DKCS +
               Religion_Muslim + Religion_Christian + Religion_Other + 
               Rural + factor(District), data=kerala, family = "binomial")


mreg6a <- glm(Owns_LPG ~ 
                Education_CWE_primary + Education_CWE_middle_high + Education_CWE_greaterhigh +
                Exp_log + NumberAdults + NumberChildrenUnder18 +
                AgeRespondent + 
                Caste_OBC + Caste_ScheduledCaste + Caste_ScheduledTribe + Caste_DKCS +
                Religion_Muslim + Religion_Christian + Religion_Other + 
                Rural + factor(District), data=kerala, family = "binomial")

mreg6 <- glm(Owns_LPG ~ Perception_Index_LPGGoodAffordable_Norm +
               Education_CWE_primary + Education_CWE_middle_high + Education_CWE_greaterhigh +
               Exp_log + NumberAdults + NumberChildrenUnder18 +
               AgeRespondent + 
               Caste_OBC + Caste_ScheduledCaste + Caste_ScheduledTribe + Caste_DKCS +
               Religion_Muslim + Religion_Christian + Religion_Other + 
               Rural + 
               factor(District), data=kerala, family = "binomial")



mreg_df <- data.frame(matrix(NA,1,5))
colnames(mreg_df) <- c("term", "estimate", "std.error", "statistic", "p.value")
mreg_df <- rbind(mreg_df,
                  tidy(mreg1a)[2:4,],
                  tidy(mreg1)[3:5,],
                  tidy(mreg2a)[2:4,],
                  tidy(mreg2)[3:5,],
                  tidy(mreg3a)[2:4,],
                  tidy(mreg3)[3:5,],
                  tidy(mreg4a)[2:4,],
                  tidy(mreg4)[3:5,],
                  tidy(mreg5a)[2:4,],
                  tidy(mreg5)[3:5,],
                  tidy(mreg6a)[2:4,],
                  tidy(mreg6)[3:5,])
mreg_df <- mreg_df[-1,]
mreg_df$model <- c("Model 1, Unadjusted", "Model 1, Unadjusted", "Model 1, Unadjusted", "Model 2, Adjusted by Perception Index", "Model 2, Adjusted by Perception Index", "Model 2, Adjusted by Perception Index")

mreg_df[mreg_df=="Perception_Index_LPGGoodAffordable_Norm"] <- "LPG is good and affordable"
mreg_df[mreg_df=="Perception_Index_FirewoodValuable_Norm"] <- "Chulha and firewood are valuable"
mreg_df[mreg_df=="Perception_Index_FirewoodUnavailableDifficult_Norm"] <- "Woodfuels are not available and inconvenient to use"
mreg_df[mreg_df=="Perception_Index_CookingHealth_Norm"] <- "Health attitudes in cooking"
mreg_df[mreg_df=="Perception_Index_GeneralHealth_Norm"] <- "General health attitudes are progressive"
mreg_df[mreg_df=="Perception_Index_ChulhaDifficult_Norm"] <- "Using the chulha is difficult"

mreg_df[mreg_df=="Education_CWE_primary"] <- "Primary School"
mreg_df[mreg_df=="Education_CWE_middle_high"] <- "Middle or High School"
mreg_df[mreg_df=="Education_CWE_greaterhigh"] <- "Greater than High School"

mreg_df$gender <- "Male"
mreg_df$estimate <- as.numeric(mreg_df$estimate)
mreg_df$std.error <- as.numeric(mreg_df$std.error)


mreg_df$perception <- c(rep("Chulha and firewood are valuable", times=6),
                         rep("Using the chulha is difficult", times=6),
                         rep("Woodfuels are unavailable, inconvenient", times=6),
                         rep("Health attitudes in cooking", times=6),
                         rep("General health attitudes are progressive", times=6),
                         rep("LPG is good and affordable", times=6))

# Figure 6: Kerala Female Ed -----

mfreg1a <- glm(Owns_LPG ~ 
                 Education_female_primary + Education_female_middle_high + Education_female_greaterhigh +
                 Exp_log + NumberAdults + NumberChildrenUnder18 +
                 AgeRespondent + 
                 Caste_OBC + Caste_ScheduledCaste + Caste_ScheduledTribe + Caste_DKCS +
                 Religion_Muslim + Religion_Christian + Religion_Other + 
                 Rural + factor(District), data=kerala, family = "binomial")

mfreg1 <- glm(Owns_LPG ~ Perception_Index_FirewoodValuable_Norm +
                Education_female_primary + Education_female_middle_high + Education_female_greaterhigh +
                Exp_log + NumberAdults + NumberChildrenUnder18 +
                AgeRespondent + 
                Caste_OBC + Caste_ScheduledCaste + Caste_ScheduledTribe + Caste_DKCS +
                Religion_Muslim + Religion_Christian + Religion_Other + 
                Rural + 
                factor(District), data=kerala, family = "binomial")


mfreg2a <- glm(Owns_LPG ~ 
                 Education_female_primary + Education_female_middle_high + Education_female_greaterhigh +
                 Exp_log + NumberAdults + NumberChildrenUnder18 +
                 AgeRespondent + 
                 Caste_OBC + Caste_ScheduledCaste + Caste_ScheduledTribe + Caste_DKCS +
                 Religion_Muslim + Religion_Christian + Religion_Other + 
                 Rural + factor(District), data=kerala, family = "binomial")

mfreg2 <- glm(Owns_LPG ~ Perception_Index_ChulhaDifficult_Norm + 
                Education_female_primary + Education_female_middle_high + Education_female_greaterhigh +
                Exp_log + NumberAdults + NumberChildrenUnder18 +
                AgeRespondent + 
                Caste_OBC + Caste_ScheduledCaste + Caste_ScheduledTribe + Caste_DKCS +
                Religion_Muslim + Religion_Christian + Religion_Other + 
                Rural + 
                factor(District), data=kerala, family = "binomial")


mfreg3a <- glm(Owns_LPG ~ 
                 Education_female_primary + Education_female_middle_high + Education_female_greaterhigh +
                 Exp_log + NumberAdults + NumberChildrenUnder18 +
                 AgeRespondent + 
                 Caste_OBC + Caste_ScheduledCaste + Caste_ScheduledTribe + Caste_DKCS +
                 Religion_Muslim + Religion_Christian + Religion_Other + 
                 Rural + factor(District), data=kerala, family = "binomial")

mfreg3 <- glm(Owns_LPG ~ Perception_Index_FirewoodUnavailableDifficult_Norm +
                Education_female_primary + Education_female_middle_high + Education_female_greaterhigh +
                Exp_log + NumberAdults + NumberChildrenUnder18 +
                AgeRespondent + 
                Caste_OBC + Caste_ScheduledCaste + Caste_ScheduledTribe + Caste_DKCS +
                Religion_Muslim + Religion_Christian + Religion_Other + 
                Rural + 
                factor(District), data=kerala, family = "binomial")

mfreg4a <- glm(Owns_LPG ~ 
                 Education_female_primary + Education_female_middle_high + Education_female_greaterhigh +
                 Exp_log + NumberAdults + NumberChildrenUnder18 +
                 AgeRespondent + 
                 Caste_OBC + Caste_ScheduledCaste + Caste_ScheduledTribe + Caste_DKCS +
                 Religion_Muslim + Religion_Christian + Religion_Other + 
                 Rural + factor(District), data=kerala, family = "binomial")

mfreg4 <- glm(Owns_LPG ~ Perception_Index_CookingHealth_Norm +
                Education_female_primary + Education_female_middle_high + Education_female_greaterhigh +
                Exp_log + NumberAdults + NumberChildrenUnder18 +
                AgeRespondent + 
                Caste_OBC + Caste_ScheduledCaste + Caste_ScheduledTribe + Caste_DKCS +
                Religion_Muslim + Religion_Christian + Religion_Other + 
                Rural + 
                factor(District), data=kerala, family = "binomial")


mfreg5a <- glm(Owns_LPG ~ 
                 Education_female_primary + Education_female_middle_high + Education_female_greaterhigh +
                 Exp_log + NumberAdults + NumberChildrenUnder18 +
                 AgeRespondent + 
                 Caste_OBC + Caste_ScheduledCaste + Caste_ScheduledTribe + Caste_DKCS +
                 Religion_Muslim + Religion_Christian + Religion_Other + 
                 Rural + 
                 factor(District), data=kerala, family = "binomial")

mfreg5 <- glm(Owns_LPG ~ Perception_Index_GeneralHealth_Norm +
                Education_female_primary + Education_female_middle_high + Education_female_greaterhigh +
                Exp_log + NumberAdults + NumberChildrenUnder18 +
                AgeRespondent + 
                Caste_OBC + Caste_ScheduledCaste + Caste_ScheduledTribe + Caste_DKCS +
                Religion_Muslim + Religion_Christian + Religion_Other + 
                Rural + factor(District), data=kerala, family = "binomial")


mfreg6a <- glm(Owns_LPG ~ 
                 Education_female_primary + Education_female_middle_high + Education_female_greaterhigh +
                 Exp_log + NumberAdults + NumberChildrenUnder18 +
                 AgeRespondent + 
                 Caste_OBC + Caste_ScheduledCaste + Caste_ScheduledTribe + Caste_DKCS +
                 Religion_Muslim + Religion_Christian + Religion_Other + 
                 Rural + factor(District), data=kerala, family = "binomial")

mfreg6 <- glm(Owns_LPG ~ Perception_Index_LPGGoodAffordable_Norm +
                Education_female_primary + Education_female_middle_high + Education_female_greaterhigh +
                Exp_log + NumberAdults + NumberChildrenUnder18 +
                AgeRespondent + 
                Caste_OBC + Caste_ScheduledCaste + Caste_ScheduledTribe + Caste_DKCS +
                Religion_Muslim + Religion_Christian + Religion_Other + 
                Rural + 
                factor(District), data=kerala, family = "binomial")


mfreg_df <- data.frame(matrix(NA,1,5))
colnames(mfreg_df) <- c("term", "estimate", "std.error", "statistic", "p.value")
mfreg_df <- rbind(mfreg_df,
                  tidy(mfreg1a)[2:4,],
                  tidy(mfreg1)[3:5,],
                  tidy(mfreg2a)[2:4,],
                  tidy(mfreg2)[3:5,],
                  tidy(mfreg3a)[2:4,],
                  tidy(mfreg3)[3:5,],
                  tidy(mfreg4a)[2:4,],
                  tidy(mfreg4)[3:5,],
                  tidy(mfreg5a)[2:4,],
                  tidy(mfreg5)[3:5,],
                  tidy(mfreg6a)[2:4,],
                  tidy(mfreg6)[3:5,])
mfreg_df <- mfreg_df[-1,]
mfreg_df$model <- c("Model 1, Unadjusted", "Model 1, Unadjusted", "Model 1, Unadjusted", "Model 2, Adjusted by Perception Index", "Model 2, Adjusted by Perception Index", "Model 2, Adjusted by Perception Index")

mfreg_df[mfreg_df=="Perception_Index_LPGGoodAffordable_Norm"] <- "LPG is good and affordable"
mfreg_df[mfreg_df=="Perception_Index_FirewoodValuable_Norm"] <- "Chulha and firewood are valuable"
mfreg_df[mfreg_df=="Perception_Index_FirewoodUnavailableDifficult_Norm"] <- "Woodfuels are not available and inconvenient to use"
mfreg_df[mfreg_df=="Perception_Index_CookingHealth_Norm"] <- "Health attitudes in cooking"
mfreg_df[mfreg_df=="Perception_Index_GeneralHealth_Norm"] <- "General health attitudes are progressive"
mfreg_df[mfreg_df=="Perception_Index_ChulhaDifficult_Norm"] <- "Using the chulha is difficult"

mfreg_df[mfreg_df=="Education_female_primary"] <- "Primary School"
mfreg_df[mfreg_df=="Education_female_middle_high"] <- "Middle or High School"
mfreg_df[mfreg_df=="Education_female_greaterhigh"] <- "Greater than High School"

mfreg_df$gender <- "Female"
mfreg_df$estimate <- as.numeric(mfreg_df$estimate)
mfreg_df$std.error <- as.numeric(mfreg_df$std.error)


mfreg_df$perception <- c(rep("Chulha and firewood are valuable", times=6),
                         rep("Using the chulha is difficult", times=6),
                         rep("Woodfuels are unavailable, inconvenient", times=6),
                         rep("Health attitudes in cooking", times=6),
                         rep("General health attitudes are progressive", times=6),
                         rep("LPG is good and affordable", times=6))

# Figure 6: Kerala combining CWE and Female ED -----

k_mediation <- mfreg_df %>%
  rbind(mreg_df)

k_mediation$p.value.bh <- p.adjust(k_mediation$p.value, method = "BH", n = length(k_mediation$p.value))
k_mediation$p.value.bh.star <- ifelse(k_mediation$p.value.bh<0.01, "***",
                                      ifelse(k_mediation$p.value.bh<0.05, "**",
                                             ifelse(k_mediation$p.value.bh<0.10, "*","")))


k_mediation_adjusted <- k_mediation %>%
  filter(model=="Model 2, Adjusted by Perception Index")

k_indices_lpg_mediation <- ggplot(k_mediation_adjusted, aes(x=factor(perception), y=estimate,
                                                    ymin=estimate - 1.96*std.error, ymax = estimate + 1.96*std.error, label=p.value.bh.star,
                                                    shape=term)) + 
  geom_point(position = position_dodge(width = .5), size=3) +
  geom_errorbar(position = position_dodge(width = .5), size=1.5, width = 0) + 
  geom_text(aes(label=p.value.bh.star, x=perception, y=estimate + 2.55*std.error), position = position_dodge(width = .5), vjust=0.75, size = 6) +
  geom_hline(yintercept=0, colour="grey60", linetype=2) + 
  theme_bw() +
  xlab("") + ylab("") +
  ggtitle("(A) Kerala: Mediation of Education-LPG Association by Perceptions") +
  coord_flip() + 
  theme(legend.position = c(-0.2,1.05),
        legend.title = element_blank(),
        legend.text = element_text(size=12),
        legend.background = element_blank(),
        legend.box.background = element_rect(colour = "black"),
        plot.title = element_text(size=16),
        axis.text = element_text(size=15),
        strip.text.x = element_text(size=14),
        strip.text.y = element_text(size=14)) +
  facet_grid(.~gender) 


k_indices_lpg_mediation

# Figure 6: Rajasthan -----


# Figure 6: Rajasthan CWE Ed -----

r_mreg1 <- glm(Owns_LPG ~ Perception_Index_FirewoodUnavailableDifficult +
                 Education_CWE_primary + Education_CWE_middle_high + Education_CWE_greaterhigh +
                 Exp_log + NumberAdults + NumberChildrenUnder18 +
                 AgeRespondent + 
                 Caste_OBC + Caste_ScheduledCaste + Caste_ScheduledTribe + Caste_DKCS +
                 Religion_Muslim + Religion_Christian + Religion_Other + 
                 Rural + 
                 factor(District), data=rajasthan, family = "binomial")


r_mreg2 <- glm(Owns_LPG ~ Perception_Index_LPGGood + 
                 Education_CWE_primary + Education_CWE_middle_high + Education_CWE_greaterhigh +
                 Exp_log + NumberAdults + NumberChildrenUnder18 +
                 AgeRespondent + 
                 Caste_OBC + Caste_ScheduledCaste + Caste_ScheduledTribe + Caste_DKCS +
                 Religion_Muslim + Religion_Christian + Religion_Other + 
                 Rural + 
                 factor(District), data=rajasthan, family = "binomial")

r_mreg3 <- glm(Owns_LPG ~ Perception_Index_LPGaffordable +
                 Education_CWE_primary + Education_CWE_middle_high + Education_CWE_greaterhigh +
                 Exp_log + NumberAdults + NumberChildrenUnder18 +
                 AgeRespondent + 
                 Caste_OBC + Caste_ScheduledCaste + Caste_ScheduledTribe + Caste_DKCS +
                 Religion_Muslim + Religion_Christian + Religion_Other + 
                 Rural + 
                 factor(District), data=rajasthan, family = "binomial")

r_mreg4 <- glm(Owns_LPG ~ Perception_Index_CookingHealth +
                 Education_CWE_primary + Education_CWE_middle_high + Education_CWE_greaterhigh +
                 Exp_log + NumberAdults + NumberChildrenUnder18 +
                 AgeRespondent + 
                 Caste_OBC + Caste_ScheduledCaste + Caste_ScheduledTribe + Caste_DKCS +
                 Religion_Muslim + Religion_Christian + Religion_Other + 
                 Rural + 
                 factor(District), data=rajasthan, family = "binomial")


r_mreg5 <- glm(Owns_LPG ~ Perception_Index_FirewoodValuable +
                 Education_CWE_primary + Education_CWE_middle_high + Education_CWE_greaterhigh +
                 Exp_log + NumberAdults + NumberChildrenUnder18 +
                 AgeRespondent + 
                 Caste_OBC + Caste_ScheduledCaste + Caste_ScheduledTribe + Caste_DKCS +
                 Religion_Muslim + Religion_Christian + Religion_Other + 
                 Rural + factor(District), data=rajasthan, family = "binomial")



# Figure 6: Rajasthan Female Ed -----


r_mfreg1 <- glm(Owns_LPG ~ Perception_Index_FirewoodUnavailableDifficult + 
                 Education_4_FemaleSpouse_primary + Education_4_FemaleSpouse_middle_high + Education_4_FemaleSpouse_greaterhigh +
                 Exp_log + NumberAdults + NumberChildrenUnder18 +
                 AgeRespondent + 
                 Caste_OBC + Caste_ScheduledCaste + Caste_ScheduledTribe + Caste_DKCS +
                 Religion_Muslim + Religion_Christian + Religion_Other + 
                 Rural + factor(District), data=rajasthan, family = "binomial")



r_mfreg2 <- glm(Owns_LPG ~ Perception_Index_LPGGood + 
                 Education_4_FemaleSpouse_primary + Education_4_FemaleSpouse_middle_high + Education_4_FemaleSpouse_greaterhigh +
                 Exp_log + NumberAdults + NumberChildrenUnder18 +
                 AgeRespondent + 
                 Caste_OBC + Caste_ScheduledCaste + Caste_ScheduledTribe + Caste_DKCS +
                 Religion_Muslim + Religion_Christian + Religion_Other + 
                 Rural + factor(District), data=rajasthan, family = "binomial")



r_mfreg3 <- glm(Owns_LPG ~ Perception_Index_LPGaffordable +
                 Education_4_FemaleSpouse_primary + Education_4_FemaleSpouse_middle_high + Education_4_FemaleSpouse_greaterhigh +
                 Exp_log + NumberAdults + NumberChildrenUnder18 +
                 AgeRespondent + 
                 Caste_OBC + Caste_ScheduledCaste + Caste_ScheduledTribe + Caste_DKCS +
                 Religion_Muslim + Religion_Christian + Religion_Other + 
                 Rural + factor(District), data=rajasthan, family = "binomial")



r_mfreg4 <- glm(Owns_LPG ~ Perception_Index_CookingHealth + 
                 Education_4_FemaleSpouse_primary + Education_4_FemaleSpouse_middle_high + Education_4_FemaleSpouse_greaterhigh +
                 Exp_log + NumberAdults + NumberChildrenUnder18 +
                 AgeRespondent + 
                 Caste_OBC + Caste_ScheduledCaste + Caste_ScheduledTribe + Caste_DKCS +
                 Religion_Muslim + Religion_Christian + Religion_Other + 
                 Rural + factor(District), data=rajasthan, family = "binomial")


r_mfreg5 <- glm(Owns_LPG ~ Perception_Index_FirewoodValuable + 
                 Education_4_FemaleSpouse_primary + Education_4_FemaleSpouse_middle_high + Education_4_FemaleSpouse_greaterhigh +
                 Exp_log + NumberAdults + NumberChildrenUnder18 +
                 AgeRespondent + 
                 Caste_OBC + Caste_ScheduledCaste + Caste_ScheduledTribe + Caste_DKCS +
                 Religion_Muslim + Religion_Christian + Religion_Other + 
                 Rural + 
                 factor(District), data=rajasthan, family = "binomial")


# Figure 6: Rajasthan combining regressions -----

mediation_regs_df <- data.frame(matrix(NA,1,5))
colnames(mediation_regs_df) <- c("term", "estimate", "std.error", "statistic", "p.value")
mediation_regs_df <- rbind(mediation_regs_df,
                  tidy(r_mreg1)[3:5,],
                  tidy(r_mreg2)[3:5,],
                  tidy(r_mreg3)[3:5,],
                  tidy(r_mreg4)[3:5,],
                  tidy(r_mreg5)[3:5,],
                  tidy(r_mfreg1)[3:5,],
                  tidy(r_mfreg2)[3:5,],
                  tidy(r_mfreg3)[3:5,],
                  tidy(r_mfreg4)[3:5,],
                  tidy(r_mfreg5)[3:5,])
mediation_regs_df <- mediation_regs_df[-1,]

mediation_regs_df$model <- c(rep("CWE",15), rep("Female Spouse",15))

mediation_regs_df$perception <- c(rep("Perception_Index_FirewoodUnavailableDifficult", 3),
                              rep("Perception_Index_LPGGood",3),
                              rep("Perception_Index_LPGaffordable",3),
                              rep("Perception_Index_CookingHealth",3),
                              rep("Perception_Index_FirewoodValuable",3))

mediation_regs_df$p.value.bh <- p.adjust(mediation_regs_df$p.value, method = "BH", n = length(mediation_regs_df$p.value))
mediation_regs_df$p.value.bh.star <- ifelse(mediation_regs_df$p.value.bh<0.01, "***",
                                             ifelse(mediation_regs_df$p.value.bh<0.05, "**",
                                                    ifelse(mediation_regs_df$p.value.bh<0.10, "*","")))


mediation_regs_df[mediation_regs_df=="Perception_Index_LPGGood"] <- "LPG is good quality cooking and beneficial"
mediation_regs_df[mediation_regs_df=="Perception_Index_LPGaffordable"] <- "LPG is desirable, but expensive"
mediation_regs_df[mediation_regs_df=="Perception_Index_CookingHealth"] <- "Cooking fuel choices affect health"
mediation_regs_df[mediation_regs_df=="Perception_Index_FirewoodUnavailableDifficult"] <- "Woodfuels are not available, inconvenient"
mediation_regs_df[mediation_regs_df=="Perception_Index_FirewoodValuable"] <- "Chulha and woodfuel are valuable"

mediation_regs_df[mediation_regs_df=="Education_CWE_primary"] <- "Primary School"
mediation_regs_df[mediation_regs_df=="Education_CWE_middle_high"] <- "Middle and High School"
mediation_regs_df[mediation_regs_df=="Education_CWE_greaterhigh"] <- "Greather than High School"

mediation_regs_df[mediation_regs_df=="Education_4_FemaleSpouse_primary"] <- "Primary School"
mediation_regs_df[mediation_regs_df=="Education_4_FemaleSpouse_middle_high"] <- "Middle and High School"
mediation_regs_df[mediation_regs_df=="Education_4_FemaleSpouse_greaterhigh"] <- "Greather than High School"


r_indices_lpg_mediation <- ggplot(mediation_regs_df, aes(x=factor(perception), y=estimate,
                                                                 ymin=estimate - 1.96*std.error, ymax = estimate + 1.96*std.error, label=p.value.bh.star,
                                                                 shape=term)) + 
  geom_point(position = position_dodge(width = .5), size=3) +
  geom_errorbar(position = position_dodge(width = .5), size=1.5, width = 0) + 
  geom_text(aes(label=p.value.bh.star, x=perception, y=((estimate + 1.96*std.error)+0.015)), position = position_dodge(width = .5), vjust=0.70, size = 6) +
  geom_hline(yintercept=0, colour="grey60", linetype=2) + 
  theme_bw() +
  xlab("") + ylab("") +
  ggtitle("(B) Rajasthan: Mediation of Education-LPG Association by Perceptions") +
  coord_flip() + 
  theme(legend.position = c(-0.2,1.05),
        legend.title = element_blank(),
        legend.text = element_text(size=12),
        legend.background = element_blank(),
        legend.box.background = element_rect(colour = "black"),
        plot.title = element_text(size=16),
        axis.text = element_text(size=15),
        strip.text.x = element_text(size=14),
        strip.text.y = element_text(size=14)) +
  facet_grid(.~model) 


r_indices_lpg_mediation


plot_grid(
  
  k_indices_lpg_mediation,
  r_indices_lpg_mediation,
  nrow=2
  
)



# Figure 7 -------

# Coefficient plot for linear regressions showing the associations between perception indices
# and years of LPG ownership in Kerala and Rajasthan households. Baseline level is \1-2 Years."
# Regressions also control for covariates and district fixed effects (not shown). Coefficients are not
# exponentiated. Whiskers show 95% confidence intervals. Stars show statistical significance after
# Benjamini-Hochberg method for controlling for false discovery rate: * P<0.10, ** P<0.05, ***
# P<0.01.

# Figure 7: Kerala -----

kerala$LPG_yearswith <- ifelse(kerala$LPG_yearswith==98, NA, kerala$LPG_yearswith)
kerala$LPG_yearswith_12 <- ifelse(kerala$LPG_yearswith==1, 1, 0)
kerala$LPG_yearswith_25 <- ifelse(kerala$LPG_yearswith==2, 1, 0)
kerala$LPG_yearswith_510 <- ifelse(kerala$LPG_yearswith==3, 1, 0)
kerala$LPG_yearswith_10plus <- ifelse(kerala$LPG_yearswith==4, 1, 0)


treg1 <- lm(Perception_Index_FirewoodValuable_Norm ~ LPG_yearswith_25 + LPG_yearswith_510 + LPG_yearswith_10plus +
               Exp_log + NumberAdults + NumberChildrenUnder18 +
               AgeRespondent + 
               Caste_OBC + Caste_ScheduledCaste + Caste_ScheduledTribe + Caste_DKCS +
               Religion_Muslim + Religion_Christian + Religion_Other + 
               Rural +
               factor(District), data=kerala)

treg2 <- lm(Perception_Index_ChulhaDifficult_Norm ~ LPG_yearswith_25 + LPG_yearswith_510 + LPG_yearswith_10plus +
               Exp_log + NumberAdults + NumberChildrenUnder18 +
               AgeRespondent + 
               Caste_OBC + Caste_ScheduledCaste + Caste_ScheduledTribe + Caste_DKCS +
               Religion_Muslim + Religion_Christian + Religion_Other + 
               Rural +
               factor(District), data=kerala)


treg3 <- lm(Perception_Index_FirewoodUnavailableDifficult_Norm ~ LPG_yearswith_25 + LPG_yearswith_510 + LPG_yearswith_10plus +
               Exp_log + NumberAdults + NumberChildrenUnder18 +
               AgeRespondent + 
               Caste_OBC + Caste_ScheduledCaste + Caste_ScheduledTribe + Caste_DKCS +
               Religion_Muslim + Religion_Christian + Religion_Other + 
               Rural + 
               factor(District), data=kerala)

treg4 <- lm(Perception_Index_CookingHealth_Norm ~ LPG_yearswith_25 + LPG_yearswith_510 + LPG_yearswith_10plus +
               Exp_log + NumberAdults + NumberChildrenUnder18 +
               AgeRespondent + 
               Caste_OBC + Caste_ScheduledCaste + Caste_ScheduledTribe + Caste_DKCS +
               Religion_Muslim + Religion_Christian + Religion_Other + 
               Rural + 
               factor(District), data=kerala)


treg5 <- lm(Perception_Index_GeneralHealth_Norm ~ LPG_yearswith_25 + LPG_yearswith_510 + LPG_yearswith_10plus +
               Exp_log + NumberAdults + NumberChildrenUnder18 +
               AgeRespondent + 
               Caste_OBC + Caste_ScheduledCaste + Caste_ScheduledTribe + Caste_DKCS +
               Religion_Muslim + Religion_Christian + Religion_Other + 
               Rural + 
               factor(District), data=kerala)

treg6 <- lm(Perception_Index_LPGGoodAffordable_Norm ~ LPG_yearswith_25 + LPG_yearswith_510 + LPG_yearswith_10plus +
               Exp_log + NumberAdults + NumberChildrenUnder18 +
               AgeRespondent + 
               Caste_OBC + Caste_ScheduledCaste + Caste_ScheduledTribe + Caste_DKCS +
               Religion_Muslim + Religion_Christian + Religion_Other + 
               Rural + 
               factor(District), data=kerala)


treg_df1 <- data.frame(matrix(NA,1,5))
colnames(treg_df1) <- c("term", "estimate", "std.error", "statistic", "p.value")
treg_df <- rbind(treg_df1,
                 tidy(treg1)[2:4,],
                 tidy(treg2)[2:4,],
                 tidy(treg3)[2:4,],
                 tidy(treg4)[2:4,],
                 tidy(treg5)[2:4,],
                 tidy(treg6)[2:4,])
treg_df <- treg_df[-1,]

treg_df[treg_df=="LPG_yearswith_25"] <- "2-5 Years"
treg_df[treg_df=="LPG_yearswith_510"] <- "5-10 Years"
treg_df[treg_df=="LPG_yearswith_10plus"] <- ">10 Years"

treg_df$perception <- c(rep("Chulha and firewood are valuable", times=3),
                        rep("Using the chulha is difficult", times=3),
                        rep("Woodfuels are unavailable, inconvenient", times=3),
                        rep("Cooking fuel choices affect health", times=3),
                        rep("General health attitudes are progressive", times=3),
                        rep("LPG is good and affordable", times=3))


years_ordered <- c("2-5 Years", "5-10 Years", ">10 Years")

treg_df$state <- "Kerala"


treg_df$p.value.bh <- p.adjust(treg_df$p.value, method = "BH", n = length(treg_df$p.value))
treg_df$p.value.bh.star <- ifelse(treg_df$p.value.bh<0.01, "***",
                                       ifelse(treg_df$p.value.bh<0.05, "**",
                                              ifelse(treg_df$p.value.bh<0.10, "*","")))


k_perception_time_fig <- ggplot(treg_df, aes(x=term, y=estimate, ymin=estimate - 1.96*std.error, ymax = estimate + 1.96*std.error)) + 
  geom_point(position = position_dodge(width = 0.3), size=3) +
  geom_errorbar(position = position_dodge(width = 0.3), size=1.5, width = 0) + 
  geom_hline(yintercept=0, colour="grey60", linetype=2) + 
  scale_x_discrete(limits=rev(years_ordered)) +
  scale_color_manual(values=c("grey50", "black")) + 
  scale_shape_manual(values=c(16,15)) +
  theme_bw() +
  xlab("") + ylab("Coefficient Estimate") +
  coord_flip() + 
  theme(axis.text.y = element_text(size=18),
        legend.text = element_text(size=14),
        axis.text.x = element_text(size=12),
        axis.title.x = element_text(size=22),
        strip.text = element_text(size=10)) +
  facet_grid(.~perception)



# Figure 7: Rajasthan  -------

rajasthan$LPG_yearswith_12 <- ifelse(rajasthan$LPG_yearswith==1, 1, 0)
rajasthan$LPG_yearswith_25 <- ifelse(rajasthan$LPG_yearswith==2, 1, 0)
rajasthan$LPG_yearswith_510 <- ifelse(rajasthan$LPG_yearswith==3, 1, 0)
rajasthan$LPG_yearswith_10plus <- ifelse(rajasthan$LPG_yearswith==4, 1, 0)
rajasthan$LPG_yearswith <- ifelse(rajasthan$LPG_yearswith==98, NA, rajasthan$LPG_yearswith)


r_perception_time_reg1 <- lm(Perception_Index_FirewoodUnavailableDifficult  ~ 
                               LPG_yearswith_25 + LPG_yearswith_510 + LPG_yearswith_10plus +
                  Exp_log + NumberAdults + NumberChildrenUnder18 +
                  AgeRespondent + 
                  Caste_OBC + Caste_ScheduledCaste + Caste_ScheduledTribe + Caste_DKCS +
                  Religion_Muslim + Religion_Christian + Religion_Other + 
                  Rural + factor(District), data = rajasthan)

r_perception_time_reg2 <- lm(Perception_Index_LPGGood  ~ 
                               LPG_yearswith_25 + LPG_yearswith_510 + LPG_yearswith_10plus +
                  Exp_log + NumberAdults + NumberChildrenUnder18 +
                  AgeRespondent + 
                  Caste_OBC + Caste_ScheduledCaste + Caste_ScheduledTribe + Caste_DKCS +
                  Religion_Muslim + Religion_Christian + Religion_Other + 
                  Rural + factor(District), data = rajasthan)

r_perception_time_reg3 <- lm(Perception_Index_LPGaffordable  ~ 
                               LPG_yearswith_25 + LPG_yearswith_510 + LPG_yearswith_10plus +
                  Exp_log + NumberAdults + NumberChildrenUnder18 +
                  AgeRespondent + 
                  Caste_OBC + Caste_ScheduledCaste + Caste_ScheduledTribe + Caste_DKCS +
                  Religion_Muslim + Religion_Christian + Religion_Other + 
                  Rural + factor(District), data = rajasthan)

r_perception_time_reg4 <- lm(Perception_Index_CookingHealth  ~ 
                               LPG_yearswith_25 + LPG_yearswith_510 + LPG_yearswith_10plus +
                 Exp_log + NumberAdults + NumberChildrenUnder18 +
                 AgeRespondent + 
                 Caste_OBC + Caste_ScheduledCaste + Caste_ScheduledTribe + Caste_DKCS +
                 Religion_Muslim + Religion_Christian + Religion_Other + 
                 Rural + factor(District), data = rajasthan)

r_perception_time_reg5 <- lm(Perception_Index_FirewoodValuable ~ 
                               LPG_yearswith_25 + LPG_yearswith_510 + LPG_yearswith_10plus +
                  Exp_log + NumberAdults + NumberChildrenUnder18 +
                  AgeRespondent + 
                  Caste_OBC + Caste_ScheduledCaste + Caste_ScheduledTribe + Caste_DKCS +
                  Religion_Muslim + Religion_Christian + Religion_Other + 
                  Rural + factor(District), data = rajasthan)


pt_reg_table <- data.frame(matrix(NA,1,5))
colnames(pt_reg_table) <- c("term", "estimate", "std.error", "statistic", "p.value")
pt_reg_table <- rbind(pt_reg_table,
                 tidy(r_perception_time_reg1)[2:4,],
                 tidy(r_perception_time_reg2)[2:4,],
                 tidy(r_perception_time_reg3)[2:4,],
                 tidy(r_perception_time_reg4)[2:4,],
                 tidy(r_perception_time_reg5)[2:4,])
pt_reg_table <- pt_reg_table[-1,]

pt_reg_table[pt_reg_table=="LPG_yearswith_25"] <- "2-5 Years"
pt_reg_table[pt_reg_table=="LPG_yearswith_510"] <- "5-10 Years"
pt_reg_table[pt_reg_table=="LPG_yearswith_10plus"] <- ">10 Years"

pt_reg_table$p.value.bh <- p.adjust(pt_reg_table$p.value, method = "BH", n = length(pt_reg_table$p.value))
pt_reg_table$p.value.bh.star <- ifelse(pt_reg_table$p.value.bh<0.01, "***",
                                       ifelse(pt_reg_table$p.value.bh<0.05, "**",
                                              ifelse(pt_reg_table$p.value.bh<0.10, "*","")))

pt_reg_table$term <- factor(pt_reg_table$term, levels = c("2-5 Years", "5-10 Years", ">10 Years"))

pt_reg_table$perception <- mediation_regs_df$perception <- c(rep("Perception_Index_FirewoodUnavailableDifficult", 3),
                                                             rep("Perception_Index_LPGGood",3),
                                                             rep("Perception_Index_LPGaffordable",3),
                                                             rep("Perception_Index_CookingHealth",3),
                                                             rep("Perception_Index_FirewoodValuable",3))
  

pt_reg_table[pt_reg_table=="Perception_Index_LPGGood"] <- "LPG is good quality cooking and beneficial"
pt_reg_table[pt_reg_table=="Perception_Index_LPGaffordable"] <- "LPG is desirable, but expensive"
pt_reg_table[pt_reg_table=="Perception_Index_CookingHealth"] <- "Cooking fuel choices affect health"
pt_reg_table[pt_reg_table=="Perception_Index_FirewoodUnavailableDifficult"] <- "Woodfuels are unavailable, inconvenient"
pt_reg_table[pt_reg_table=="Perception_Index_FirewoodValuable"] <- "Chulha and firewood are valuable"

pt_reg_table$state <- "Rajasthan"


perception_time_table_combined <- pt_reg_table %>%
  bind_rows(treg_df)

perception_time_table_combined$term <- factor(perception_time_table_combined$term, levels = c("2-5 Years","5-10 Years", ">10 Years"))


ggplot(perception_time_table_combined, aes(x=perception, y=estimate,
                         ymin=estimate - 1.96*std.error, ymax = estimate + 1.96*std.error, label=p.value.bh.star,
                         shape=term)) + 
  geom_point(position = position_dodge(width = 0.5), size=3) +
  geom_errorbar(position = position_dodge(width = 0.5), size=1.5, width = 0) + 
  geom_text(aes(label=p.value.bh.star, x = perception, y=((estimate-1.96*std.error)-0.015)), 
            position = position_dodge(width = .5),  vjust=0.7, size = 6) +
  geom_hline(yintercept=0, colour="grey60", linetype=2) + 
  theme_bw() +
  xlab("") + ylab("") +
  coord_flip() + 
  theme(legend.title = element_blank(),
        legend.text = element_text(size=14),
        legend.position = c(0.8,0.94),
        axis.text = element_text(size=18),
        strip.text.x = element_text(size=10),
        strip.text.y = element_text(size=14)) + 
  facet_wrap(state~.)







